#coding=utf-8
import math
import re
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from utils.mylib_utils import *

# 原位气球：
# AL_TKQQ_TPWX_2023080823.txt
def evaluate_uncertainty_balloon(L2_filepath, system_param, k=1, given_param="U"):
    print("评定参数1：", system_param)
    print("评定参数2：", given_param)

    L2_filename = os.path.basename(L2_filepath)
    info = L2_filename.split("_")
    # given_param = "T"               # 传参（测试）
    device_code = info[1]           # 设备编码：探空气球等
    data_type = info[2]
    data_L2 = parse_QQ_data(L2_filepath)    # 高度(距地)   气温   气压   风向   风速
    #print(data_L2)
    #if device_code == "TKQQ" and data_type == "TPWX":
    if data_type == "TPWX":
        if given_param == "T":
            data_L2 = data_L2[:, (0, 1, 4)]
        elif given_param == "P":
            data_L2 = data_L2[:, (0, 2)]
        elif given_param == "W":
            data_L2 = data_L2[:, (0, 4)]
        # elif given_param == "X":
            # data_L2 = data_L2[:, (0, 3)]

        param = given_param

    #print(data_L2)
    #return
    result_df = None
    # try:
    if param == "T":
        result_df = get_GUM_Balloon_T_uncertainty(
            data_L2,
            k,                                     # GUM参数k
            system_param["T_error"],               # 温度载荷的最大允许误差，表示正负0.3（输入）
            system_param["k_square"],              # 系统所在海拔高度, km
            system_param["param_u2"],              # 用户输入参数（输入）
            system_param["Vd"],                    # 探空仪上升的速度（输入）（m/s）
            system_param["cp"],                    # 空气的定压比热（单位：kJ/(kg·K)）
            system_param["t"],                     # 传感器响应时间（输入）（s）
            system_param["t_high_ratio"],          # 温度随高度的变化率（输入）
            system_param["param_u5"],              # 探空仪传感结构引入的测量不确定度分量
            system_param["t0"]                     # 探空仪数据采样周期（s）
        )
    elif param == "P":
        result_df = get_GUM_Balloon_P_uncertainty(
            data_L2,
            k,                                     # GUM参数k
            # xxxxx
        )
    elif param == "W":
        result_df = get_GUM_Balloon_W_uncertainty(
            data_L2,
            k,                                     # GUM参数k
            system_param['deltaV']                 # 风速系数
        )
    # except:
    #     print('File Error: '+L2_filepath)
    return result_df


# 温度不确定度计算
def get_GUM_Balloon_T_uncertainty(
        data_L2,
        K,
        T_error = 0.3,      # 温度载荷的最大允许误差，表示正负0.3（输入）
        k_square = 3,       # sqrt_k计算系数
        param_u2 = 0,       # 用户输入参数（输入）
        Vd = 8,             # 探空仪上升的速度（输入）（m/s）
        cp = 1.004,         # 空气的定压比热（单位：kJ/(kg·K)）
        t = 3,              # 传感器响应时间（输入）（s）
        t_high_ratio = "(0,11]km:-0.007;(11,20]km:0.00;(20,32]km:0.001",  #温度随高度的变化率（输入）
        param_u5 = 0,
        t0 = 1,             # 探空仪数据采样周期（s）
        # s = 30,             # s为高度间隔（m）（内部计算获得）
):
    # 参数：
    # height = 30        # 输入：修改为向量
    #
    # T_error = 0.3      # 温度载荷的最大允许误差，表示正负0.3（输入）
    # k_square = 3       # u1计算系数
    # u2 = 0             # 用户输入参数（输入）
    # Vd = 8             # 探空仪上升的速度（输入）（m/s）
    # V = 8              # 该高度所测风速（输入）
    # cp = 1.004         # 空气的定压比热（单位：kJ/(kg·K)）
    # t = 3              # 传感器响应时间（输入）（s）
    # t_high_ratio = "(0,11]km:-0.007;(11,20]km:0.00;(20-32]km:0.001"  #温度随高度的变化率（输入）
    # u5 = 0
    # t0 = 1             # 探空仪数据采样周期（s）
    # s = 30             # s为高度间隔（m）（内部计算获得）

    # 获取高度列，温度列
    z = data_L2[:, 0] / 1000        # 高度单位转换为 km
    T = data_L2[:, 1]
    W = data_L2[:, 2]
    shape = np.shape(z)
    # ratio = np.full(np.shape(H), np.nan)

    # 计算中使用的常量
    heigh_diff = (z[1] - z[0]) * 1000    # 按照米计算高度差
    k_sqrt = math.sqrt(k_square)

    # 不确定度分量
    # 1.由温度传感器自身准确度引入的测量不确定度分量计算
    u1 = np.full(shape, T_error / k_sqrt)

    # 2.由太阳辐射引入的测量不确定度分量计算
    u2 = np.full(shape, param_u2)

    # 3.由气动效应引入的测量不确定度分量计算（W为对应高度的风速）
    u3 = ((Vd**2 + W**2) / (2 * k_sqrt * cp)) * 0.001

    # 4.由传感器响应时间引入的测量不确定度分量计算
    t_high_ratio_arr = parseParamStr(t_high_ratio)
    dtdz = getParamByHeightV2(z, t_high_ratio_arr)
    u4 = np.abs(((t * Vd) / k_sqrt) * dtdz)
    u4 = np.abs(u4)

    # 5.由探空仪传感结构引入的测量不确定度分量计算
    u5 = np.full(shape, param_u5)
    u5 = np.abs(u5)

    # 6.由数据传输引入的测量不确定度分量计算
    u6 = ((t0 * Vd) / k_sqrt) * dtdz
    u6 = np.abs(u6)

    # 7.由数据处理引入的测量不确定度分量计算
    u7 = (heigh_diff / k_sqrt) * dtdz
    u7 = np.abs(u7)

    # 8.计算合成标准不确定度
    # 合成标准不确定度
    u = np.sqrt(
        u1**2
        + u2**2
        + u3**2
        + u4**2
        + u5**2
        + u6**2
        + u7**2
    )
    # 扩展不确定度
    U = K * u
    # 不确定度
    U, u, u1, u2, u3, u4, u5, u6, u7 = np.round(
        (U, u, u1, u2, u3, u4, u5, u6, u7), 4
    )
    result = np.vstack((z, T, U, u, u1, u2, u3, u4, u5, u6, u7)).T
    header = [
        "高度(km)",
        "温度(K)",
        "U(K)_GUM扩展不确定度",
        "u(K)_标准不确定度",
        "u1(K)_由温度传感器自身准确度引入",
        "u2(K)_由太阳辐射引入",
        "u5(K)_由气动效应引入",
        "u6(K)_由传感器响应时间引入",
        "u8(K)_由探空仪传感结构引入",
        "u9(K)_由数据传输引入",
        "u10(K)_由数据处理引入",
    ]
    result_df = pd.DataFrame(result, columns=header)
    return result_df


# 压力不确定度计算
def get_GUM_Balloon_P_uncertainty(
        data_L2,
        K,
        k_square = 3,       # sqrt_k计算系数
        Vd = 8,             # 探空仪上升的速度（输入）（m/s）
        t = 3,              # 传感器响应时间（输入）（s）
        t_high_ratio = "(0,4]km:-12;(4,11]km:-8;(11,20]km:-4;(20,32]km:-1",  #气压随高度的变化率（输入）
        t0 = 1,             # 探空仪数据采样周期（s）
):
    # 获取高度列，压力
    z = data_L2[:, 0] / 1000        # 高度单位转换为 km
    P = data_L2[:, 1]
    shape = np.shape(z)

    # 计算中使用的常量
    heigh_diff = (z[1] - z[0]) * 1000    # 按照米计算高度差
    k_sqrt = math.sqrt(k_square)

    # 不确定度分量
    # 1.由压力传感器自身准确度引入的测量不确定度分量计算
    deltaP = getSensorErrorValueByPressure(P)
    u1 = deltaP / k_sqrt

    # 2.由传感器温度特性引入的测量不确定度分量计算
    u2 = np.full(shape, 0)

    # 3.由传感器气动效应引入的测量不确定度分量计算
    u3 = np.full(shape, 0)

    # 4.由传感器响应时间引入的测量不确定度分量计算
    t_high_ratio_arr = parseParamStr(t_high_ratio)
    dpdz = getParamByHeightV2(z, t_high_ratio_arr)
    u4 = np.abs(((t * Vd) / k_sqrt) * dpdz)

    # 5.由数据传输的测量不确定度分量计算
    u5 = ((t0 * Vd) / k_sqrt) * dpdz
    u5 = np.abs(u5)

    # 6.由数据处理引入的测量不确定度分量计算
    u6 = (heigh_diff / k_sqrt) * dpdz
    u6 = np.abs(u6)

    # [计算合成标准不确定度]
    # 合成标准不确定度
    u = np.sqrt(
        u1**2
        + u2**2
        + u3**2
        + u4**2
        + u5**2
        + u6**2
    )
    # [扩展不确定度]
    U = K * u
    # 返回
    U, u, u1, u2, u3, u4, u5, u6 = np.round(
        (U, u, u1, u2, u3, u4, u5, u6), 4
    )
    result = np.vstack((z, P, U, u, u1, u2, u3, u4, u5, u6)).T
    header = [
        "高度(km)",
        "气压(Pa)",
        "U(K)_GUM扩展不确定度",
        "u(K)_标准不确定度",
        "u3(K)_由压力传感器自身准确度引入",
        "u4(K)_由温度特性引入",
        "u5(K)_由气动效应引入",
        "u6(K)_由传感器响应时间引入",
        "u9(K)_由数据传输引入",
        "u10(K)_由数据处理引入",
    ]
    result_df = pd.DataFrame(result, columns=header)
    return result_df


# 风速不确定度计算
def get_GUM_Balloon_W_uncertainty(
        data_L2,
        K,
        k_square = 3,       # u1计算系数
        deltaV = 2          # 风速系数，用户输入，默认2m/s
):
    # 参数：
    # height = 30        # 输入：修改为向量
    #
    # T_error = 0.3      # 温度载荷的最大允许误差，表示正负0.3（输入）
    # k_square = 3       # u1计算系数
    # u2 = 0             # 用户输入参数（输入）
    # Vd = 8             # 探空仪上升的速度（输入）（m/s）
    # V = 8              # 该高度所测风速（输入）
    # cp = 1.004         # 空气的定压比热（单位：kJ/(kg·K)）
    # t = 3              # 传感器响应时间（输入）（s）
    # t_high_ratio = "(0,11]km:-0.007;(11,20]km:0.00;(20-32]km:0.001"  #温度随高度的变化率（输入）
    # u5 = 0
    # t0 = 1             # 探空仪数据采样周期（s）
    # s = 30             # s为高度间隔（m）（内部计算获得）

    # 获取高度列，温度列
    z = data_L2[:, 0] / 1000        # 高度单位转换为 km
    W = data_L2[:, 1]
    shape = np.shape(z)
    # ratio = np.full(np.shape(H), np.nan)

    # 计算中使用的常量
    # heigh_diff = (z[1] - z[0]) * 1000    # 按照米计算高度差
    k_sqrt = math.sqrt(k_square)

    # 不确定度分量
    # 1.由气球速度测量引入的测量不确定度分量计算
    u1 = 0.001 * W

    # 2.由大气风速解算引入的测量不确定度分量计算
    deltaP = getSensorErrorValueByWindSpeed(W)
    u2 = deltaP / k_sqrt

    # 3.由气球风阻系数引入的测量不确定度分量计算
    u3_tmp = deltaV / k_sqrt
    u3 = np.full(shape, u3_tmp)

    # [计算合成标准不确定度]
    # 合成标准不确定度
    u = np.sqrt(
        u1**2
        + u2**2
        + u3**2
    )
    # [扩展不确定度]
    U = K * u
    # 精度处理
    U, u, u1, u2, u3 = np.round(
        (U, u, u1, u2, u3), 4
    )
    result = np.vstack((z, W, U, u, u1, u2, u3)).T
    header = [
        "高度(km)",
        "温度(K)",
        "U(K)_GUM扩展不确定度",
        "u(K)_标准不确定度",
        "u11(K)_由气球速度测量引入",
        "u12(K)_由大气风速解算引入",
        "u13(K)_由气球风阻系数引入",
    ]
    result_df = pd.DataFrame(result, columns=header)
    return result_df


# 根据风速判断参数
def getSensorErrorValueByWindSpeed(W):
    t_high_ratio_arr = [('-', 10, 1), (10, '-', 0.1)]
    ratio = np.full(np.shape(W), np.nan)
    for ratio_param in t_high_ratio_arr:
        if ratio_param[0] == '-':
            mask = (W < ratio_param[1])
            ratio[mask] = ratio_param[2]
        elif ratio_param[1] == '-':
            mask = (ratio_param[0] <= W)
            ratio[mask] = ratio_param[2] * W[mask]
            #ratio[mask] = ratio_param[2]
    return ratio


# 根据气压值判断压力传感器的最大误差
# 气压≥500hPa，△(P)≤2hPa
# 气压＜500hPa，△(P)≤1hPa
#
def getSensorErrorValueByPressure(P):
    t_high_ratio_arr = [('-', 500, 1), (500, '-', 2)]
    ratio = np.full(np.shape(P), np.nan)
    for ratio_param in t_high_ratio_arr:
        if ratio_param[0] == '-':
            mask = (P < ratio_param[1])
            ratio[mask] = ratio_param[2]
        elif ratio_param[1] == '-':
            mask = (ratio_param[0] <= P)
            ratio[mask] = ratio_param[2]
    return ratio


# 根据高度判断参数值
# 参数格式：
# t_high_ratio_arr = [(0, 11, -0.007), (11, 20, 0.00), (20, 32, 0.001)]
# t_high_ratio_arr可使用函数parseParamStr获取
def getParamByHeightV2(H, t_high_ratio_arr):
    ratio = np.full(np.shape(H), np.nan)
    for ratio_param in t_high_ratio_arr:
        mask = (ratio_param[0] <= H) & (H <= ratio_param[1])
        ratio[mask] = ratio_param[2]
    return ratio


# def getParamByHeightV3(z):
#     t_high_ratio = [(0, 11, -0.007), (11, 20, 0.00), (20, 32, 0.001)]
#
#     height_km = z / 1000
#     # height_km = np.array(z) / 1000
#     # height_km = z.apply(lambda x: x / 1000)
#     ratios = np.array(t_high_ratio)
#     # conditions = (height_km > ratios[:, 0]) & (height_km < ratios[:, 1])
#     # conditions = [((x > ratios[:, 0]) & (x < ratios[:, 1])).any() for x in zip(*height_km)]
#     conditions = np.any((height_km[:, None] > ratios[:, 0]) & (height_km[:, None] < ratios[:, 1]), axis=1)
#     # values = np.where(conditions, ratios[:, 2], 0)
#     # values = np.where(conditions, np.repeat(ratios[:, 2], len(height_km)), 0)
#     values = np.where(conditions, np.repeat(ratios[:, 2], ratios[:, 1] - ratios[:, 0]), 0)
#     return values.max()

##########################
# 解析分段函数参数为数组
# 输入："(0,11]km:-0.007;(11,20]km:0.00;(20,32]km:0.001"
# 输出：[(0, 11, -0.007), (11, 20, 0.00), (20, 32, 0.001)]
def parseParamStr(piecewise_fun_str):
    ret = []
    arr = re.split(";", piecewise_fun_str)
    for item in arr:
        arr1 = re.split(":", item)
        if len(arr1) >= 2:
            value = arr1[1]         # value = -0.007
            print(value)
            # 使用正则表达式进行匹配和提取
            flag = False
            pattern = r'\((.*?)\]'  # 匹配以括号开头，中间任意字符，以]结尾的部分。支持： (a, b]
            match = re.search(pattern, arr1[0])
            if match:
                flag = True
                condition_text = match.group(1)  # 获取匹配结果中的第一个捕获组（codition_text: 0,11; 11,20; 20,32）
                # print(condition_text)
            else:
                pattern = r'\((.*?)\)'  # 支持：(a, b)
                match = re.search(pattern, arr1[0])
                if match:
                    flag = True
                    condition_text = match.group(1)
            if flag:
                arr3 = re.split(",", condition_text)
                print(arr3)
                element = (float(arr3[0]), float(arr3[1]), float(value))
                ret.append(element)
            else:
                print("分段函数参数错误（level-2）！");
                continue
        else:
            print("分段函数错误（level-1）！")
            continue
    return ret


# 未使用（废弃）
# 解析参数
# "(0,11]km:-0.007;(11,20]km:0.00;(20,32]km:0.001"
def getParamByHeight(height, param_str):
    height_km = float(height)/1000
    arr = re.split(";", param_str)
    for item in arr:
        arr1 = re.split(":", item)
        if len(arr1) >= 2:
            value = arr1[1]         # value = -0.007
            print(value)
            # 使用正则表达式进行匹配和提取
            pattern = r'\((.*?)\]'  # 匹配以括号开头，中间任意字符，以]结尾的部分
            match = re.search(pattern, arr1[0])
            if match:
                condition_text = match.group(1)  # 获取匹配结果中的第一个捕获组（codition_text: 0,11; 11,20; 20,32）
                # print(condition_text)
                arr3 = re.split(",", condition_text)
                print(arr3)
                if len(arr3)<2:
                    print("分段函数参数错误!");
                    continue
                if height_km > float(arr3[0]) and height_km <= float(arr3[1]):
                    return value
            else:
                print("No match found.")
                continue
        else:
            print("分段函数参数错误！")
            continue
    return 0


def testGetParamByHeight():
    param1 =  "(0,11]km:-0.007;(11,20]km:0.00;(20,32]km:0.001"
    height = "22000"
    ret = getParamByHeight(height, param1)
    print(ret)
def testParseParamStr():
    param1 =  "(0,11]km:-0.007;(11,20]km:0.00;(20,32]km:0.001"
    ret = parseParamStr(param1)
    print(ret)


def parse_QQ_data(filepath):
    if not os.path.exists(filepath):
        print('Filename Error: {}'.format(filepath))
        return 0
    head_index = get_head_index_QQ(filepath)
    filename = os.path.basename(filepath)
    fileinfo = filename.split('_')
    device_code = fileinfo[1]       # 设备编码：探空气球等
    data_type = fileinfo[2]

    data = np.empty(0)
    #if device_code == "TKQQ" and data_type == "TPWX":
    if data_type == "TPWX":
        encode = get_file_encod(filepath)
        probe_data = pd.read_csv(filepath, skiprows=head_index, encoding=encode, sep='\s+')
        data = probe_data.values
        # 默认情况：
        #      0           1       2        3          4        5       6        7           8        9
        #   高度(距地)   时间    气温      气压      湿度     露点    温露差    虚温       风向     风速
        #      0         0.0     25.9       857        5     -16.6     42.5    299.3       221        2
        # 存在其他情况：
        #      0           1       2     3    4     5        6
        #   高度(距地)    气温   气压  湿度  风向  风速     虚温
        if "时间" not in probe_data.head():
            data = data[:, (0, 1, 2, 4, 5)]
        else:
             data = data[:, (0, 2, 3, 8, 9)]
    else:
        print('Filename Error: {}'.format(filepath))

    # plt.figure()
    # plt.plot(data[:, 1:], data[:, 0])
    # plt.ylabel('Height/km')
    # plt.show()
    data[data == -9999] = np.nan
    return data


def get_head_index_QQ(filepath):
    encode = get_file_encod(filepath)
    with open(filepath, 'r', encoding=encode) as f:
        filehead_all = f.readlines()
    head_index = 0
    while not (filehead_all[head_index].strip().startswith('Rang') or filehead_all[head_index].strip().startswith('Alt') or
               filehead_all[head_index].strip().startswith('#End of Header') or filehead_all[head_index].strip().startswith('高度(距地)')):
        head_index += 1
    return head_index


def get_file_encod(file_path):
    # 检查文件编码
    encode = get_file_encoding(file_path)
    print("file_encode:", encode)
    if encode == "utf-8":
        read_file_encoding = 'utf8'
    else:
        read_file_encoding = 'GBK'
    return read_file_encoding


if __name__ == "__main__":
    L2_filepath = r"C:\tmp\test_data_From_YJ_2311\原位探测new\探空气球数据-5个文件\AL_TKQQ_TPWX_2023080823.txt"
    # L0_path = r"D:\1.Work\1.Data\新型地基\数据\实测数据\20230920\2023年09月20日数据分享\L0"
    ret = evaluate_uncertainty_balloon(L2_filepath, {})
    print("======================")
    print(ret)

    # 测试
    # testGetParamByHeight()
    # testParseParamStr()